
* British Election Study
/*
Use the following variables:
cutsTooFarNational
privatTooFar
howToReduceDeficit
enviroProtection
immigrationLevel
redistSelf
euRefVote
govtHandouts

*/



clear all

capture program drop store_est
program define store_est
	local indepvar = "`1'"
	matrix mat_`indepvar'_n[`2',`3']=e(N)
	matrix mat_`indepvar'_beta[`2',`3']=_b[temp_x]
	matrix mat_`indepvar'_se[`2',`3']=_se[temp_x]
	matrix mat_`indepvar'_var[`2',`3']=_se[temp_x]^2
	matrix mat_`indepvar'_t[`2',`3']=_b[temp_x]/_se[temp_x]
end

* SE of effect estimate average
capture program drop calc_SE
program define calc_SE

	local indepvar = "`1'"
	local model = `2'

	* Total sum of variances
	egen sum_var_`indepvar'_`model'=sum(mat_`indepvar'_var`model')
	
	* Second term
	local tempsum = 0
	local actualn = 0
	forvalues i=1(1)`3' {
		if mat_`indepvar'_se[`i',`model']!=. {
			local actualn = `actualn'+1
		}
		forvalues j=1(1)`3' {
			if `i'!=`j' & mat_`indepvar'_se[`i',`model']!=. & mat_`indepvar'_se[`j',`model']!=. {
				local tempsum = `tempsum' + `4'corr_mat[`i',`j'] * mat_`indepvar'_se[`i',`model'] * ///
					mat_`indepvar'_se[`j',`model'] * `4'flip_`indepvar'[`i']  * `4'flip_`indepvar'[`j']
			}
		}
	}
	
	* Total standard error
	gen se_`indepvar'_`model' = sqrt(`tempsum' + sum_var_`indepvar'_`model'[1])*(1/`actualn') in 1

end

cd "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Quantifying bias\replication_code"
use "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Quantifying bias\ext_data\bes_panel_W4.dta" 

* Dependents
order cutsTooFarNational privatTooFar howToReduceDeficit enviroProtection immigrationLevel redistSelf euRefVote govtHandouts
replace euRefVote=. if euRefVote==2 

local j=1
foreach i of varlist cutsTooFarNational - govtHandouts  {
	replace `i' = . if `i' == 9999 | `i' == 99
	gen outcome`j' = `i'
	label variable outcome`j' "`i'"
	local j=`j'+1
	drop `i'
}


* Gen independents
* 15 or under
gen educationYears=.
	replace educationYears=9 if p_education_age==1
	replace educationYears=10 if p_education_age==2
	replace educationYears=11 if p_education_age==3
	replace educationYears=13 if p_education_age==4
	replace educationYears=15 if p_education_age==5
	replace educationYears=17 if p_edlevel==5
gen bin_edu=.
	replace bin_edu=0 if educationYears<15
	replace bin_edu=1 if educationYears>14 & educationYears!=.
gen income=.
	replace income = 2500 if p_gross_personal == 1
	replace income = 7500 if p_gross_personal == 2
	replace income = 12500 if p_gross_personal == 3
	replace income = 17500 if p_gross_personal == 4
	replace income = 22500 if p_gross_personal == 5
	replace income = 27500 if p_gross_personal == 6
	replace income = 32500 if p_gross_personal == 7
	replace income = 37500 if p_gross_personal == 8
	replace income = 42500 if p_gross_personal == 9
	replace income = 47500 if p_gross_personal == 10
	replace income = 55000 if p_gross_personal == 11
	replace income = 65000 if p_gross_personal == 12
	replace income = 85000 if p_gross_personal == 13
	replace income = 115000 if p_gross_personal == 14
	replace income = income/1000

* Gen controls
gen sex = gender
gen municipality = pcon
gen occupation = ns_sec // kanske ns_sec_analytic


* Inclusion
gen include_indicator=.
replace include_indicator=1 if sex!=. & age!=. & municipality!=. & occupation!=. & income!= . & educationYears != .

* Cleanup
keep outcome1 - include_indicator age


* Create necessary matrices for storage
foreach indepvar in bin_edu educationYears income {
	matrix mat_`indepvar'_beta=J(8,2,.)
	matrix mat_`indepvar'_n=J(8,2,.)
	matrix mat_`indepvar'_se=J(8,2,.)
	matrix mat_`indepvar'_var=J(8,2,.)
	matrix mat_`indepvar'_t=J(8,2,.)
}

* Correlation matrix
pwcorr outcome1-outcome8 if include_indicator==1
matrix corr_mat=r(C)


* Run analyses
local i=1
foreach depvar of varlist outcome1-outcome8 {
	
	foreach indepvar in educationYears bin_edu income {
	
		di `i' ", `depvar' `indepvar'"
	
		qui:su `depvar' if `indepvar'!=. & include_indicator!=.
		qui:gen temp_y = (`depvar'-r(mean))/r(sd)
		qui:su `indepvar' if `depvar'!=. & include_indicator!=.
		qui:gen temp_x = (`indepvar'-r(mean))/r(sd)
	
		* Empty model
		qui: reghdfe temp_y temp_x c.sex##i.age if include_indicator==1, noabsorb
		store_est `indepvar' `i' 1
	
		* Naive model
		if "`indepvar'" == "bin_edu"|"`indepvar'" == "educationYears" {
			qui: reghdfe temp_y temp_x c.sex##i.age i.occupation income if include_indicator==1, absorb(municipality)
			store_est `indepvar' `i' 2
		}
		else {
			qui: reghdfe temp_y temp_x c.sex##i.age i.occupation educationYears if include_indicator==1, absorb(municipality)
			store_est `indepvar' `i' 2
		}
		
		drop temp_x temp_y
		
	}
	local i=`i'+1
}

drop _all
foreach indepvar in bin_edu educationYears income {

	* Turn matrices into variables
	svmat mat_`indepvar'_n
	svmat mat_`indepvar'_beta
	svmat mat_`indepvar'_t
	svmat mat_`indepvar'_se
	svmat mat_`indepvar'_var
	
	* Flip vector
	gen flip_`indepvar' = 1
	replace flip_`indepvar' = -1 if mat_`indepvar'_beta1<0
	
	* Flipped betas
	gen `indepvar'_1 = mat_`indepvar'_beta1 * flip_`indepvar'
	gen `indepvar'_2 = mat_`indepvar'_beta2 * flip_`indepvar'
	
	* Means
	egen m_`indepvar'_1 = mean(`indepvar'_1)
	egen m_`indepvar'_2 = mean(`indepvar'_2)
	gen `indepvar'_m =.
	replace `indepvar'_m = m_`indepvar'_1[1] in 1
	replace `indepvar'_m = m_`indepvar'_2[1] in 2
	drop m_*
	
	* Standard errors
	calc_SE `indepvar' 1 8
	calc_SE `indepvar' 2 8
	gen `indepvar'_se =.
	replace `indepvar'_se = se_`indepvar'_1[1] in 1
	replace `indepvar'_se = se_`indepvar'_2[1] in 2

	gen `indepvar'_Nmin = .
	gen `indepvar'_Nmax = .
	gen `indepvar'_delta = .
	gen `indepvar'_tdelta = .
	forvalues i = 1(1)2 {
		egen Nmin_`indepvar'_`i' = min(mat_`indepvar'_n`i')
		egen Nmax_`indepvar'_`i' = max(mat_`indepvar'_n`i')
		replace `indepvar'_Nmin = Nmin_`indepvar'_`i'[1] in `i'
		replace `indepvar'_Nmax = Nmax_`indepvar'_`i'[1] in `i'
		drop Nmin_`indepvar'_`i' Nmax_`indepvar'_`i'
	}
	replace `indepvar'_delta = `indepvar'_m[2] - `indepvar'_m[1] in 2
	replace `indepvar'_tdelta = `indepvar'_delta[2]/sqrt(`indepvar'_se[1]^2 + `indepvar'_se[2]^2) in 2
	
}

keep *_m *_se *_Nmin *Nmax *_delta *_tdelta
keep if [_n]<3
rename * BES_*
save BES_est.dta, replace


